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ABSTRACT 

The general consensus is that in order to reproduce the observed solar p-mode oscillation frequencies, 
turbulence should be included in solar models. However, until now there has not been any well-tested 
efficient method to incorporate turbulence into solar modeling. We present here two methods to include 
turbulence in solar modeling within the framework of the mixing length theory, using the turbulent 
velocity obtained from numerical simulations of the highly superadiabatic layer of the sun at three 
stages of its evolution. The first approach is to include the turbulent pressure alone, and the second 
is to include both the turbulent pressure and the turbulent kinetic energy. The latter is achieved by 
introducing two variables: the turbulent kinetic energy per unit mass, and the effective ratio of specific 
heats due to the turbulent perturbation. These are treated as additions to the standard thermodynamic 
coordinates (e.g. pressure and temperature). We investigate the effects of both treatments of turbulence 
on the structure variables, the adiabatic sound speed, the structure of the highly superadiabatic layer, 
and the p-mode frequencies. We find that the second method reproduces the SAL structure obtained in 
3D simulations, and produces a p-mode frequency correction an order of magnitude better than the first 
method. 

Subject headings: turbulence — Sun: interior 



1. INTRODUCTION 

The effects of turbulence on the structure of a solar 
model depends not only on how it is modeled, but also on 
how it is incorporated into the solar model. The change 
in the p-mode oscillation frequencies caused by turbulence 
can be used as a measure of the effects of turbulence on 
the structure of solar models, since helioseismology pro- 
vides an opportunity to probe directly and sensitively solar 
structure. In the standard solar models (SSM), the tur- 
bulent convection is modeled in terms of the local mixing- 
length theory (MLT), but the turbulent pressure, turbu- 
lent kinetic energy and turbulent entropy are ignored. The 
fact that the computed frequencies of p-modes from stan- 
dard solar models are higher than the observed values, 
shows that standard solar models need to be refined. The 
frequency dependence of the discrepancy reveals that there 
is a problem with the model, and that the problem lies 
very near the surface. Balmforth (1992a, 1992b, 1992c) 
uses the non-local mixing-length theory to model the tur- 
bulent convection, and includes the turbulent pressure, 
but ignores the turbulent kinetic energy and turbulent en- 
tropy in modeling the sun. The computed frequencies of p- 
modes from such a model are even higher than those com- 
puted from the standard solar models (Balmforth 1992a). 
Canuto (1990,1996) has proposed a semi-analytical model 
of the turbulent convection. The main idea is to include a 
full turbulent spectrum in the model of convection. Appli- 
cations to stellar models have been discussed by Canuto & 
Mazzitelli (1991,1992) and Canuto, Goldman, & Mazzitelli 
(1996). The parameters of the Canuto-Mazzitelli (C&M) 
approach are derived from the laboratory experiments of 



the incompressible convection, and extrapolated to the 
stellar conditions. Using this approach, the superadiabatic 
peak is much higher than that of the standard solar mod- 
els, while the computed frequencies of solar p-modes are 
closer to the observed values than those from the standard 
solar models. 

Turbulence is a highly nonlinear phenomenon. Numer- 
ical experiments are a direct, effective way to investi- 
gate turbulence, in addition to the laboratory experiments 
(Niemela et al 2000). Chan & Sofia (1989) performed a 
three-dimensional numerical simulation of the deep com- 
pressible convection, where the superadiabatic excess in 
the temperature gradient is very small. They showed that 
in the deep regions, the mixing- length approximation is 
valid. They derived expressions for key physical parame- 
ters such as the convective energy flux. These, in principle, 
can be conveniently applied to stellar models. The Chan- 
Sofia formula for the convective flux was incorporated into 
the Yale stellar evolution code by Lydon (1993) and Ly- 
don, Fox, & Sofia (1992). The diffusion approximation was 
used to treat radiative transfer, and the radiative atmo- 
sphere was treated by Lydon et al. (1992) as in the stan- 
dard MLT Yale models (e.g., Guenther et al. 1992). They 
found that the peak of the superadiabatic layer (SAL) is 
not as high as in the C&M model and is located at a greater 
depth below the surface. However, their models yielded a 
larger discrepancy between the observed solar p-mode fre- 
quencies and those computed from standard solar models. 
Part of the discrepancy could be attributed to the fact 
that these models did not match the solar radius very pre- 
cisely. But this is expected since the simulation, which is 
valid for the deeper adiabatic regions, is extrapolated into 
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the SAL, where the temperature gradient greatly exceeds 
the adiabatic gradient, and some inaccuracy is inevitable. 

In order to overcome these difficulties, Kim (1993) and 
Kim ct al. (1995, 1996) conducted a three-dimensional nu- 
merical simulation whose domain includes shallower layers. 
This simulation treats the coupling of radiation and con- 
vection and includes realistic equation of state and radia- 
tive opacities, taken from the Yale Stellar Evolution Code. 
The Kim et al.'s (1995, 1996) simulation treats radiative 
transfer in the diffusion approximation, which is only valid 
in an optically thick medium, and consequently cannot be 
used in the solar atmosphere. Thus the top boundary of 
the simulation was set below the SAL peak. The Kim et 
al.'s (1996) models were parameterized as a varying mixing 
length with depth by Demarque, Guenther, & Kim (1997), 
in precisely calibrated solar models. The p-mode frequen- 
cies of the models were found to agree more closely with 
the observed solar p-mode frequencies, than the standard 
solar model frequencies, but these models exhibit, like the 
C&M models, a higher SAL peak than the standard solar 
models. 

More recently, Kim & Chan (1998) have completed a 
three-dimensional radiative hydrodynamic simulation of 
the complete extent of the SAL including the solar at- 
mosphere (about 2 pressure scale heights above and 2.5 
pressure scale heights below the SAL). The numerical ap- 
proach was described by Kim & Chan (1997). The simu- 
lation of Kim & Chan is fully compressible and uses the 
realistic equation of state and opacities. The radiation has 
been treated by utilizing the three-dimensional Edding- 
ton approximation, which is valid in both the optically 
thin regions near the surface and the optically thick re- 
gions in the interior. Demarque, Guenther, & Kim (1999) 
mimicked this simulation in the calibrated solar models 
by increasing the opacity coefficient k near the surface. 
Such models are called perturbed-K models. The discrep- 
ancy between calculated and observed p-mode frequencies 
decreases when compared to standard solar models. A 
deeper layer (about 5 pressure scale heights above and 6 
pressure scale heights below the SAL) was simulated by 
Stein & Nordlund (1998) using different numerical meth- 
ods for the convective and radiative components. Instead 
of parameterizing the simulation, Rosenthal et al. (1999) 
matched the simulation to an envelope which was con- 
structed using an MLT envelope code. The computed fre- 
quencies are in better agreement with the observed solar 
p-mode frequencies when compared to the standard solar 
models. Abbett et al. (1997) have also discussed the same 
transition layer. 

One way to include turbulence in stellar models is to use 
the numerical simulations to directly calculate the convec- 
tive temperature gradient and its derivatives needed in 
stellar model calculations. However, because turbulence is 
chaotic, nonlocal, and three-dimensional, and because it 
involves nonlinear interactions over many disparate length 
scales, the simulation for the whole convective zone is too 
computationally expensive for stellar model calculations. 
Besides, the accuracy of the simulation, if performed, is too 
low to match the required accuracy for the stellar model 
calculations. Fortunately, the existing simulations show 
that the MLT prediction deviates from simulations only 
in the SAL part of the convection zone. Even so, it is still 



impractical to calculate the convective temperature gradi- 
ent and its derivatives, from numerical simulation of the 
SAL in the stellar evolution model caculations. Robinson 
et al. (2000) have performed some hydrodynamical simu- 
latons (the viscosity parameter c M = 0.2y/2, see below) of 
the SAL at three stages in the solar evolution: the zero- 
age main sequence (ZAMS), the present sun and the sub- 
giant. The results show that the turbulent velocities in the 
SAL, as functions of gas pressure, change little for all these 
stages. Therefore, it is feasible to compute the effects of 
turbulence on the convective gradients. Overshooting was 
observed in both the solar subatmosphere and 3D numer- 
ical simulations, and the K-modcls (Demarque et al 1999) 
demonstrate that overshooting is likely to be one of the 
keys to match the real Sun. It is a challenge to include 
overshooting within the framework of the local MLT since 
overshooting is a nonlocal phenomenon. 

In this paper, we use the results of the 3D numeri- 
cal simulations of turbulence to calculate the convective 
temperature gradient, and its derivatives needed in solar 
model calculations. As MLT is valid in most of the con- 
vective zone, it is convenient to include turbulence within 
the framework of MLT. In §2 we summarize the new re- 
sults ((^ = 0.2) of the 3D simulations of the solar SAL, 
at 3 stages of its evolution. These simulations are similar 
to those by Robinson et al. (2000), but with lower vis- 
cosities. §3 describes how to calculate the convective tem- 
perature gradient using turbulent velocities. We describe 
the calibrated solar models with turbulence in §4. In §5 
the influence of turbulence on the structure variables, the 
adiabatic sound speed, the structure of the highly supera- 
diabatic layer, and the p-mode frequencies are calculated 
and compared to the observed solar p-mode frequencies. 
Concluding remarks are presented in the last section. 



2. TURBULENT VELOCITIES 

We incorporate the radiative hydrodynamical simula- 
tions of the outer layers of the sun into the ID stellar 
models. Three 3D simulations have been performed us- 
ing the hydrostatic ID stellar models, at three stages of 
its evolution (ZAMS, present sun and subgiant), as start- 
ing points. The physics (thermodynamics, the equation 
of state, and opacities) in the simulation is the same as in 
the ID stellar models. These simulations follow closely the 
approach described by Kim & Chan (1998), and are de- 
scribed in more detail by Robinson et al. (2000). The full 
hydrodynamical equations were solved in a thin subsection 
of the stellar model, i.e. a 3D box located in the vicinity 
of the photosphere. For the radiative transport, the diffu- 
sion approximation was used in the deep region (r > 10 3 ) 
of the simulation, while the 3D Eddington approximation 
was used (Unno & Spiegel 1966) in the region above. Af- 
ter the simulation had reached a steady state, statistical 
integrations were performed for each simulation for over 
2500 seconds in the case of the solar surface convection. 

Turbulence can be measured by the turbulent Mach 
number M. = v"/v s , where v" is the turbulent velocity, 
and v s is the sound speed. The MLT is valid when M is 
sufficiently small. In the outer layers of the sun M. can be 
of order unity (Cox and Giuli 1968), but in the deep con- 
vection region M is almost zero. The turbulent velocity 
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is denned by the velocity variance: 

v'i = (^-vfy\ (i) 

where the overbar denotes a combined horizontal and tem- 
poral average, and Vi is the total velocity. Fig. 1 shows the 
run of M. as a function of log P (in base 10) in the con- 
vection simulations for the ZAMS model, the present sun 
and the subgiant model, respectively. We note that the 
maximum of M. is about 0.7 and changes little from the 
ZAMS to the present sun. Using A4, we can define the 
turbulent kinetic energy per unit mass x as 

X =\M 2 v\. (2) 
The turbulent contribution to the entropy is 

Sturb - X/T, (3) 
where T is the gas temperature. 
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Fig. 1. — Turbulent Mach number as a function of depth at 3 
evolutionary stages. 
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Fig. 2. — Ratio of turbulent to total pressure in the outer layers 
as a function of depth at 3 evolutionary stages. 

Turbulence in the stratified layers of the solar convec- 
tion zone is not isotropic. For convenience, we define the 
parameter 7 to reflect the anisotropy of turbulence, 

Pturb = (7 - 1)PX, (4) 

where p\ is the turbulent kinetic energy density. Since 
-Pturb = pv" 2 ', we can relate 7 to the turbulent velocity as 



//\2 



7=1 + 2(v'Uv") 



(5) 



7 = 5/3 when turbulence is isotropic (v" = v" = v'y); 
7 = 3 or 7 = 1 when turbulence is completely anisotropic 
(v" — v" or v" — 0, respectively). The physical meaning of 
7 is the specific heat ratio due to turbulence. This affects 
the radial turbulent pressure distribution. Fig. 2 shows 
-Pturb, in which the radial turbulent pressure is scaled with 
the gas pressure, P gas - The total pressure is defined as 



rad 



turb- 



(6) 



Note that P t urb/-P g a S = «M) 2 since v s = (P gas / p) 1/2 ■ 
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Fig. 3. — Turbulent kinetic energy per unit mass as a function 
of depth at 3 evolutionary stages. 
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Fig. 4. — Specific heat ratio due to turbulence as a function of 
depth at 3 evolutionary stages. 

The turbulent contribution to the pressure, kinetic en- 
ergy and entropy, can all be expressed in terms of x an d 
7. Therefore, turbulence can be parameterized by these 
two parameters. Figs. 3 and 4 show their variations as 
functions of depth. 
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3. CONVECTIVE TEMPERATURE GRADIENTS WITH THE 
TURBULENT VELOCITIES 

Abbett et al. (1997) tried to include turbulence in solar 
modeling within the framework of MLT by using simulated 
pressure and density in calculating the temperature gradi- 
ent in the convection zone. As they pointed out, such an 
application of MLT is not self-consistent. Lydon and Sofia 
(1995) developed a self-consistent method to include mag- 
netic fields in calculating the convective temperature gra- 
dients within the framework of MLT. Then, Lydon, Guen- 
ther and Sofia (1996) used it to successfully explain the 
observed variation of solar p-modes with the solar cycle. 
Recently, Li and Sofia (2001) have updated this method to 
reproduce the observed cyclic variations of all solar global 
parameters such as solar luminosity, solar effective tem- 
perature and solar radius. We briefly show how the same 
method can be used to calculate the influence of turbulence 
on the temperature gradients in the convection zone. 

3.1. Turbulent variables 

Introducing x, and 7, computed from the 3D simula- 
tions, the equation of state becomes 

p = p(P T ,T,x,7). (7) 

Including the turbulent kinetic energy, the first law of ther- 
modynamics becomes 

dQ T = dU + PdV + d X 

= dU T + [Pt - (7 - l)(x/V)]dV, (8) 

where Ut = U + x, is the total internal energy per unit 
mass and V = 1/p is the volume per unit mass. 

3.2. Convective stability criterion 

The convective stability criterion can still be expressed 
by the difference between density derivatives of a mass 
element and its surroundings: 

(dp/dr) e - (dp/dr) s > 0. (9) 

However, since we are using Pt, T, \ an d 7 as the inde- 
pendent variables, we have 

dp/p = pdP T /P T - p'dT/T - vdxlx - v'dlh, (10) 

where 



( d\np \ 
\ d In P T J 



I'' 



( dlnp \ 
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Pt,T,j V 17111 '/ P T ,T, X 

As a result, the stability criterion becomes 

V ra d<Vad-An. (11) 

In this expression, A m reflects the direct influence of tur- 
bulence, defined by 

A m = (i/V x + i/V 7 )V ad /M, (12) 

where we have assumed (dx/dr) e = (dx/dr) s and 
(dj/dr) e = (d-f/dr) s , while 



^ d In P T J , 



V7 , — 3 nL r P T yj _ P T 8 
V rad — i 67mcG M r T 4 V ad — pTc p 

The other symbols have the usual meanings. 
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3.3. Convective temperature gradients 

3.3.1. Flux conservation with turbulence 

The convective temperature gradient V con v is deter- 
mined by the requirement that the total energy flux -Ftotai 
equals the sum of the radiative flux F ra d and the convective 
flux F 

liuA A convi 



F 



total 



F 



rad 



F,, 



(13) 



The total flux at any given layer in the star is determined 
by the photon luminosity L r , 

T- AacG T A M r _ 

(14) 



-Ftotai 



Aitr 2 



3 KPrr 2 



rad- 



The radiative flux is determined by the convective tem- 
perature gradient: 



-Frad 



AacG T 4 M r 



V c 



(15) 



3 nP T r 2 

The convective flux is determined by the convective veloc- 
ity fconv and the heat excess DQt'- 



F„ 



PWconv DQ T - 



(16) 



When the convective velocity is much smaller than the 
sound speed of the medium, the process can be considered 
to be of constant pressure. In this case, the heat excess 
can be obtained from the first law given by Eq. (8): 



DQ T = c p DT 



P T p'v 



1 



D X 



P T p!v' 
PM7 



D 1 . (17) 



3.3.2. Mixing length approximation 

Using the mixing length approximation, DT, Dx and 
D-f can be expressed by the mixing length l m as follows: 

DT/T = (l/T)d(DT)/dr(l m /2) 

= (Vconv - V e )(W2)(Vflp), 

D X /X = (l/ X )d(D X )dr(l m /2) = 0, (18) 
D7/7 - (l/ 7 )d(D 7 )dr(l m /2) = 0, 

where H p = —PT^dr / 'cZPt) is the pressure scale height. In 
order to determine w CO nv, the MLT assumes that half of 
the work done by half of the radial buoyancy force acted 
over half the mixing length goes into the kinetic energy of 
the element (i^ onv /2). Since the radial buoyancy force per 
unit mass is related to the density difference by: 

k r - -g{Dp/p), (19) 

and since the process is in pressure equilibrium, we obtain 

l 2 m 9 



= p' (V C onv - V e ) 



8H, D ' 



(20) 



where g = GM r /r 2 is the gravitational acceleration. 
An additional relation is required to close the MLT: 

(dQT/dr) e — (radiative losses) + {change of x}, (21) 

which can be expressed as: 

{2acT 3 )/(pc p v conv )(Lu/(l + i^ 2 )(V conv - V e ) 

= (V e - V ad ) + (V ad /A*)(^V x + l/V 7 ), (22) 

where lo = npl m . 
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3.3.3. Result 
Solving Eqs. (13) and (22), we obtain 

V conv = V ad + {y/Vj 2 C)(l + y/V) - A m , (23) 
where y is the solution of the following equation 

2Ay 3 + Vy 2 + V 2 y - V = 0. (24) 
A, 7q, C, and V are defined by 

A=(9/8)[^ 2 /(3 + c 2 )], 
7o = [( Cp p)/(2a C r 3 )][(l + ^ 2 )H, 
C = (g/l 2 m »')/8H p , 

V = l/[ 7 oC 1/2 (V rad - V ad + A m f/ 2 ]. 

From these formulas it can be seen that the effect of turbu- 
lence on the temperature gradient in the convection zone 
can be taken into account by modifying the adiabatic gra- 



dient: 



V ad = [1 - (^V x + z/V^/^Vad. 



(25) 



Therefore, it is easy to include this effect in the standard 
stellar structure codes. 

4. CALIBRATED SOLAR MODELS 

4.1. Standard solar model 

For the purpose of comparison, we construct a standard 
solar model with the Yale Stellar Evolution Code. The 
OPAL opacities tables (Iglesias & Rogers 1996) are used 
together with the low-temperature opacities from Kurucz 
(1991). The equation of state is taken from Rogers, Swen- 
son, & Iglesias (1996). When out of the table, the Yale 
standard implementation with the Debye-Hukel correction 
is used (Guenther ct al.1992). Helium and heavy clement 
diffusion processes are included in the model. Heavy ele- 
ment diffusion is implemented by assuming that all heavy 
elements diffuse with the same velocity as fully ionized iron 
(Guenther & Demarque 1997). The model atmosphere 
is constructed using the empirical Krishna-Swamy (KS) 
T — t relation (Guenther et al. 1992). The solar model is 
evolved from the zero-age main sequence to the current so- 
lar age. The mixing-length ratio a and the helium content 
Y have been adjusted in the usual way so as to match the 
solar luminosity, radius, and the ratio of heavy elements 
to hydrogen {Z/X = 0.0230) at the surface of the model 
(Grevesse & Sauval 1998). These are obtained by choos- 
ing (Z , X ) and a to be (0.0188, 0.7091724) and 2.135772. 
The standard (or reference) solar model is abbreviated as 
the SSM. 

4.2. Solar model with turbulent pressure alone 

The simplest way to taken into account turbulence in so- 
lar modeling is to include turbulent pressure (or Reynolds 
stress) alone, as done by many authors (e.g., Balmforth 
1992a). In this case, only the hydrostatic equilibrium 
equation needs to be modified as follows (method 1): 

dP GM r ,„ ^ 

9mT-w (1 + /?) ' (26) 

where P = P gas + P-ad, and 

2-Pturb dP tur b\ ( - 9Pt 



pgr 



dP 



1 + 



turb 



dp 



Here 2P tur b / {pgr) originates from the spherical coordinate 
system adopted, representing a kind of geometric effect. 
The equations that govern the envelope integrations also 
need to be changed accordingly. 

We implement this case by modifying the Yale Stel- 
lar Evolution Code and obtain a nonstandard model in 
the same way we obtain the standard model. We as- 
sume Pturb, set equal to its value for the present sun, 
does not change from the ZAMS to the present age of 
the sun. The adjustable parameters now are fixed as 
(Z ,X ,a)=(0.0188, 0.7092889, 2.138190). This model is 
called the Pressure Solar Model (PSM). 

4.3. Solar model with \ and 7 as independent variables 

The form of the continuity equation and the equation of 
transport of energy by radiation is not affected by turbu- 
lence. The hydrostatic equation includes a Reynolds stress 
term due to turbulence 

dP GM r 1 d , o 

^ = -— P-^(r^), (28) 

where P = P gas + P ra d- Since the last term can be rewrit- 
ten as dPturb/dr + 2(7 — l)x/r, this equation becomes 

dP T GM r 2(7 -l)x 



dM r 



Anr A 



47rr 3 



(29) 



The last term on the right hand side of Eq. (29) also em- 
bodies the same spheric geometric effect as 2P tur b/ {pgr) 
in Eq. (27). The energy conservation equation is also af- 
fected by turbulence because the first law of thermody- 
namics should include the turbulent kinetic energy 



dL r _ _ dS T 
dM r ~ € dt ' 



(30) 



(32) 



(27) 



where 

TdS T = c p dT -^dP T +(l + ^) d X + ^dj. 

P V PPX J PPH 

(31) 

The equation of energy transport by convection, 
dT T GM r 
dW r ~ ~P7 4^ conv ' 
does not change in form, but the convection temperature 
gradient, obtained in the previous section, is different from 
that without turbulence. The equations that govern en- 
velope integrations also need to be changed accordingly. 
This method will be referred below as method 2. 

We implement this case by modifying the Yale Stellar 
Evolution Code. Once again, we obtain a nonstandard 
model in the same way as we obtained the standard model. 
If we assume that \ an d 7 do not change with time (letting 
them equal their values at the present age of the sun), the 
adjustable parameters now must be set as (Zq,Xq,(x) = 
(0.0188,0.7092715,2.271540). We use the spline interpo- 
lation of x and 7 given in §2 for their pressure dependence 
in the model calculations. We call this model the Energy 
Solar Model 1 (ESM1). 

In order to investigate the evolutionary effects of \ an d 
7, we linearly interpolate between the two simulations that 
are closest to the evolutionary state of the model. In this 
case, {Z ,X ,a) = (0.0188,0.7092945,2.271462) in order 
to match the observed L Q , Rq, and {Z/X). We call this 
model as Energy Solar Model 2 (ESM2). 
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5. INFLUENCE OF TURBULENCE ON THE SOLAR 
STRUCTURE 

We shall now investigate how different methods for in- 
cluding turbulence in solar modeling, affect the solar model 
structure. 

5.1. Measured by structural variables 

Fig. 5 depicts turbulence-induced relative changes (with 
respect to the SSM) of the solar structural variables for the 
PSM (dotted line), ESM1 (dashed line) and ESM2 (solid 
line overlapped on the dashed line), as functions of the log- 
arithm of the total pressure with base 10. The changes for 
the pressure, temperature, and luminosity are calculated 
at the same radius coordinate, while the radius change is 
calculated at the same interior mass coordinate M r . From 
this figure it can be seen that 

1. In all the cases turbulence affects the distribution 
of the pressure, temperature, and the (radiative 
and convective, rather than the total) luminosity 
around the SAL peak (specified by the dashed 
vertical line in the figure); 

2. Both methods produce almost the same maximum 
pressure change of 15%, but method 1 produces 

a larger change for the temperature (8% vs 4%), 
radiative and convective luminosity (130% vs 70% 
and 100% vs 80%); 

3. It is surprising that turbulence near the surface 
affects (slightly) the mass distribution (denoted by 
the radius change) in the core of the solar models. 

The increase of the radiative luminosity does not neces- 
sarily imply the increase of the convective gradient since 
the radiative flux depends on not only the convective gra- 
dient, but also the temperature T, the radiative opacity 
k, and the pressure Pt, as expressed by Eq. 15. Simi- 
larly, the decrease of the convective luminosity does not 
necessarily imply the increase of the convective gradient 
since the convective flux depends on the specific heat at 
constant pressure c p , density p, temperature T, convective 
velocity w CO nv, and the mixing length l m : 



-fconv — "conv 
yP 'm 



(33) 



Fig. 6 shows how p, k, v conv , l m , and c p in the PSM and 
ESM change with respect to the SSM. The increase of the 
radiative luminosity below the SAL peak is mostly gener- 
ated by the decrease of the radiative opacity caused by the 
decrease of the temperature and density. The decrease of 
the convective luminosity above the SAL peak is mostly 
generated by thedecrease of the convective velocity. 

Figs. 5 and 6 show that the temperature, density and 
total pressure decrease when the turbulent pressure is in- 
cluded. The decrease of the total pressure is possible since 
what supports gravity in the solar interior is the total pres- 
sure gradient, not the pressure itself. Besides, the gas 
pressure may decrease in order to maintain hydrodynamic 
equilibrium when turbulence provides a turbulent pres- 
sure. Consequently, the gas density and/or the gas tem- 
perature should decrease. However, the two methods make 
difference here: the decrease of the temperature (density) 



in the ESM is smaller (larger) than that in the PSM. The 
cause is that the mixing length in the ESM increases near 
the SAL peak, but the mixing length in the PSM decreases 
near the SAL peak. This implies that the transport of 
energy by convection near the peak in the ESM is more 
efficient than that in the PSM. As a result, different SAL 
structures are expected, as addressed in the next subsec- 
tion. 




Fig. 5. — Turbulence-induced relative changes of the solar struc- 
tural variables (pressure, temperature, radiative luminosity, convec- 
tive luminosity, and radius) for the PSM and ESM with respect to 
the SSM. The vertical line indicates the location of the SAL peak 
of the SSM: LogP T = 5.09. 

The comprehensive effect of turbulence on solar struc- 
ture manifests itself in the adiabatic sound speed 



AlnC= i(Alnr! + AlnP T - Amp), (34) 



where fi is the first adiabatic exponent defined by 



<91nP T 
<91np 



(35) 



S,xn 



when turbulence is modeled by \ an d 7- Fig. 7 shows 
these four quantities for the ESM and PSM. Obviously, 
the pressure and density changes contribute to the change 
of the adiabatic sound speed, as does the change of the first 
adiabatic exponent. This shows that the thermodynamic 
properties of the solar matter changes when turbulence is 
present, as expected. Nevertheless, the two methods pro- 
duce a difference for the first adiabatic exponent around 
the SAL peak. 
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Fig. 6. — Realtive change of the density, opacity, convective ve- 
locity, mixing length, and specific heat at the constant pressure in 
the PSM (dotted), and ESM (solid) with respect to the SSM. 
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Fig. 7. — Relative change of the adiabatic sound speed and its 
contributors in the ESM and PSM. 

Another feature that we can see from Figs. 5-7 [the 
relative change for any variable X is defined as follows: 
— A ssm )/A ssm ] is that although turbulence is 
restricted to the highly superadiabatic layer (log Ft € 
[4.6,7], see Figs. 1-4), its influence extends deeply into 



the solar interior for the ESMs. For example, we still 
see some influence near the base of the convective zone 
at log Pt ~ 13. This is a natural consequence of continu- 
ity. 

5.2. Measured by superadiabaticity 

The superadiabaticity (V — V a d) as a function of loga- 
rithm of total pressure with base 10 is depicted in Fig. 8 
for the SSM (solid line), and ESM1 (dotted line). The 
PSM has the same SAL as the SSM, and ESM2 has the 
same SAL as the ESM1. The SAL peak of the SSM equals 
to 0.45, while that of the ESM equals to 0.40, about 11% 
lower than that of the SSM. 

The corresponding 3D simulations produce an SAL very 
similar to the ID model. The maximum of V — V a d is 
about 0.4. The 3D solar surface simulations by Stein and 
Nordlund (Rosenthal et al. 1999) and Canuto's ID turbu- 
lence models found a much higher peak of about 0.8. The 
simulations by Stein and Nordlund employ a hyperviscos- 
ity model for the subgrid scales, while our 3D simulations 
use the Smagorinsky model (Smagorinsky 1963) with the 
smallest viscosity that was numerically stable (c M = 0.2, 
see below). If we increased the viscosity by an order of 
magnitude then the SAL peak was similar to Stein and 
Nordlund's. 




Fig. 8. — Structure of the highly superadiabatic layer (SAL). The 
SAL of the SSM is overlapped on that of the PSM. 

In order to understand these results, we note that the 
actual temperature gradient is determined by the relative 
efficiency between the radiative and convective transport 
of heat since the total energy flux is fixed by the total 
luminosity, see Eq. (14). Fig. 9 shows the actual and adi- 
abatic temperature gradients for the ESM and PSM, from 
which we can see that the SAL peak decrease in the ESM 
with respect to the SSM is mostly caused by the decrease 
of the actual gradient due to the inclusion of the turbulent 
kinetic energy. The maximum relative change of the con- 
vective (and adiabatic) gradient in the ESM is less than 
10%. In the MLT approach, the convective flux is pro- 
portional to the convective velocity, as shown by Eq. (16). 
In 3D simulations, the total flux is equal to the radiative 
flux 

frad: plus the enthalpy flux F e , plus the turbulent ki- 
netic energy flux F^. Near the SAL peak, Fk « 0. The 
enthalpy flux (Chan and Sofia 1989) is proportional to the 



8 



root mean square fluctuation of vertical velocity v" defined 
in Eq. (1), 

F e = pc~^T = pc p C[v' z TyX', (36) 

where C[v' z T'\ is the correlation coefficient between the 
temperature and vertical velocity fluctuations, and T" is 
the root mean square fluctuation of temperature. Since 
Fk ~ near the SAL peak, the decrease of the convec- 
tive gradient in the ESM and the 3D simulations implies 
an increase of the convective velocity w C onv and the rms 
fluctuation of the vertical velocity just near the peak, as 
confirmed by Fig. 10. This shows that the SAL peak is 
sensitive to the turbulent velocity. 



0.8 
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Fig. 9. — Convective and adiabatic temperature gradients. Those 
of the SSM and PSM are overlapped on each other. 

The most reliable method to determine the turbulent ve- 
locity is by direct numerical simulations (DNS) using the 
real solar kinematic viscosity v (Landau & Lifshitz 1987). 
Since the number of degrees of freedom needed to repre- 
sent 3D turbulent convection is proportional to Re 9 / 4 , to 
resolve numerically all the scales in the solar convection 
zone (where the Reynolds number Re w 10 12 ) would re- 
quire 10 27 grid points (Canuto 2000). However, the max- 
imum number of grid points allowed by the present tech- 
nology in DNS is o(10 9 ). This forces us to use large eddy 
simulations (LES) by increasing the kinematic viscosity v 
so that it represents the effects of Reynolds stresses on 
the unresolved or sub-grid scales (SGS). In the LES sim- 
ulations, we used the SGS formula due to Smagorinsky 
(1963), 

v={c^f{2a:a) 1 ' 2 . (37) 

The colon inside the brackets denots tensor contraction of 
the rate of strain = (X7 i Vj + \7jV i )/2, A is the grid spac- 
ing, and c M is an adjustable dimcnsionlcss parameter. In 
the Smagorinsky model, v w (c^A) 2 ?//^, where the char- 
acteristic u and length L arc estimated from the resolved 
motions. Consequently, the Reynolds number reads 

Re = uL/v = (L/c M A) 2 . (38) 

In nondimensional units, L ss 1, c M = 0.2 and A ss 0.01. 
So Re ~ 10 6 . This number is still much smaller than 
the Reynolds number 10 12 in the solar convection zone. 
Therefore, the viscosity must be overestimated, although 
this Reynolds number should be high enough for the fluid 



to be turbulent. The fact that the standard MLT gener- 
ates a lower convective velocity may imply that the MLT 
assumes a larger viscosity than the 3D simulations. This 
consideration disfavors those solar models with a higher 
SAL peak than that of the SSM. 
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Fig. 10. — Radial turbulent velocity and convective velocity as a 
function of depth (measured by logarithm of pressure with base 10). 
Those of the SSM and PSM are overlapped on each other. 

It should be pointed out that we compare different solar 
models at the same radius coordinate while we compare 
the SAL at the same pressure coordinate. Therefore, we 
observe a large difference for the convective velocity in 
Figs. 6 but a small (or no) difference between the ESM 
(PSM) and the SSM in Fig. 9. 

5.3. Measured by p-mode oscillation frequencies 

Our principal goal is to investigate how the treatment of 
turbulence in solar modeling affects the computed model 
structure. Therefore, we do not include the contribution 
of turbulence to the pulsation equations when we calculate 
p-mode oscillation frequency differences caused by turbu- 
lence. We use Guenther's pulsation code (1994) to calcu- 
late the p-mode frequencies under the adiabatic approxi- 
mation, for our SSM, PSM, ESM1 and ESM2, respectively. 

Fig. 11 shows the frequency differences (turbulent so- 
lar model, including the PSM, ESM1 and ESM2, minus 
standard solar model) scaled by the mode mass Q n \ (e.g., 
Christensen-Dalsgaard et al. 1991). The frequency dif- 
ferences between the PSM and SSM are comparable with 
Balmforth's (Table 1, 1992a). Both models use turbu- 
lent pressure alone. The frequency diffcrencics between 
the EMS and SSM are comparable with Rosenthal et al.'s 
(Figure 5, 1999). 

In order to examine if the frequency differences shown in 
Fig. 11 are caused by the calibration of the solar models, 
we calculated the corresponding uncalibrated models for 
the calibrated PSM and ESM2, respectively. Unlike the 
calibrated model, the starting model for the uncalibrated 
model is the standard solar model at the present age of 
the sun, instead of the ZAMS model. We switch on the 
turbulence and evolve the model 10 timesteps with 1 year 
as the timestep length. The frequency differences between 
calibrated and uncalibrated turbulent models arc less than 
1 ^Hz. 
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Fig. 11. — P-mode frequency difference diagrams. Turbulent 
model minus standard model, for the turbulent pressure solar model 
(psm), and the solar model with the turbulent pressure and turbu- 
lent kinetic energy (esml and esm2). Plotted are the 1 = 0, 1, 2, 3, 
4, 10, 20, 100 p-modes. 

By calculating the frequency differences between the 
ESM1 and ESM2 models we can find out the effect of tem- 
poral change of turbulence in the evolutionary timescale 
of the sun. From Fig. 1 1 we can see that this effect is very 
small (less than 0.5 /zHz). 



x 
=1 



15 



10 - 



Adiabatic 



psm: dot-dashed 

ssm: dashed 
esm 1 : dotted 
esm2: solid 




-5 

1500 2000 2500 3000 3500 4000 4500 
z/ e [>Hz] 

Fig. 12. — P-mode frequency difference diagrams, observation 
minus model, scaled by the mode mass Q n i, for the standard solar 
model (ssm), the turbulent pressure solar model (psm), the solar 
model with fixed turbulent pressure and kinetic energy (esml), and 
the solar model with evolutionary turbulent pressure and kinetic 
energy (esm2, almost overlap with esml). Plotted are the 1 = 0, 1, 
2, 3, 4, 10, 20,...,100 p-modcs. 

As pointed out ealier, the PSM model is obtained by 
including turbulent pressure alone in the solar modeling, 
while the ESMs are obtained by introducing the turbu- 
lent variables x an d 7 to include both turbulent pressure 
and kinetic energy. Therefore, the frequency differences 
between these two kinds of models reflects the different 
treatments of turbulence in the solar modeling. Physically, 
the differences are caused by the addition of turbulent ki- 
netic energy. Fig. 11 shows that the frequency differences 
caused by turbulent kinetic energy are much larger than 
those caused by turbulent pressure alone in size. Fig. 12 
shows that the frequency changes caused by turbulent ki- 



netic energy make the computed model frequencies match 
the solar data better than the SSM model, which is in 
agreement with Rosenthal et al.'s result (see their Fig- 
ures 1 and 6, 1999). 

6. CONCLUDING REMARKS 

We have shown how different treatments of turbulence 
in solar modeling affect the model structure within the 
framework of the standard mixing length theory. The 
turbulent velocity is obtained from 3D numerical simu- 
lations of turbulence in the highly superadiabatic layer 
of the sun. When we introduce the turbulent kinetic en- 
ergy per unit mass \ an d the effective ratio of specific 
heats due to the turbulent perturbation 7, as independent 
thermodynamic variables, the resultant solar model is in 
agreement with the patched solar model, in which the sim- 
ulated SAL replaces the original SAL. The frequency shift 
tends to match the observations better than the standard 
solar model. In contrast, when we use only the turbulent 
pressure, the turbulent effects are substantially underesti- 
mated (in the sense that the resultant p-mode frequency 
shift is much smaller). 

Another difference between method 1 and 2 is that the 
SAL peak in the ESM is lower than that of the SSM, but 
that of the PSM is the same as that of the SSM. The rea- 
son is that the increase of the mixing length, in the vicinity 
of the SAL peak, by the turbulent kinetic energy, is more 
than double the reduction of the mixing length by the tur- 
bulent pressure (see panel 4 of Fig. 6) . The SAL produced 
by method 2 is consistent with that of our 3D simulations, 
while method I does not change the SAL structure of the 
SSM. 

Previously, the turbulent pressure was considered to 
play an important role in solar modeling, but we have 
shown that this is not true: it is the turbulent kinetic en- 
ergy that is important. In fact, if we calibrate the solar 
model, the elevation caused by the turbulent pressure Ar 
vashishes. Consequently, 



Ar/ c ph 
f^dr/C 



(39) 



vanishes, where c p h is the adiabatic sound speed at the 
solar surface (photosphere). However, we obtain almost 
the same frequency changes with and without calibrating 
the model radius to the solar radius at the present age 
of the sun. Moreover, the solar model with the turbulent 
pressure alone does not reproduce the p-model frequency 
shift obtained by Rosenthal et al.'s patched model (1999). 

To match the observed solar p-mode frequencies is not 
the only motivation to improve the solar model by includ- 
ing turbulence. More interesting is to generate a solar 
model that can excite the observed p-modes and damp 
the unobserved modes simultaneously. Such an effort is in 
progress (Li et al 2001). Overshooting may also play an 
important role in this regard. 
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